Unbiased Reprogramming with 48 Factors by Reprogram-Seq

Jialei Duan

2022-02-01

---

Duan, J., Li, B., Bhakta, M., Xie, S., Zhou, P., Munshi, N.V., and Hon, G.C. (2019). Rational Reprogramming of Cellular States by Combinatorial Perturbation. Cell Rep. 27, 3486–3499.e6.

• BioProject Accession: PRJNA483190
• GEO Accession: GSE117795
• GitHub repository: jlduan/Reprogram-Seq@GitHub

---

Rendered at 2022-02-01 13:16:29 America/Chicago.

Load required packages.

library(tidyverse)
library(Matrix)
library(patchwork)
library(extrafont)

Sys.time()

[1] "2022-02-01 13:16:30 CST"

Sys.timezone()

[1] "America/Chicago"

Data preparation

Functions loading

source(
    file = file.path(
        SCRIPT_DIR,
        "utilities.R"
    )
)

Data loading

PROJECT_DIR <- file.path(
    "/Users/jialei/Maestral/Data/Projects/UTSW/Cellular_reprogramming",
    "Cardiac_reprogramming/Data_submission/Github"
)

gene_symbols <- vroom::vroom(
    file = file.path(
        PROJECT_DIR, "data", "misc", "genes.tsv"
    ),
    col_names = FALSE
)

gene_symbols <- setNames(
    object = gene_symbols$X2,
    nm = gene_symbols$X1
)

gene_symbols |> head()

ENSMUSG00000051951 ENSMUSG00000089699 ENSMUSG00000102343 ENSMUSG00000025900 
            "Xkr4"           "Gm1992"          "Gm37381"              "Rp1" 
ENSMUSG00000109048 ENSMUSG00000025902 
             "Rp1"            "Sox17" 

length(gene_symbols)

[1] 27999

Matrix

matrix_readcount_use <- Matrix::sparseMatrix(
    i = readRDS(
        file.path(
            PROJECT_DIR, "data/drop-seq", "expr_readcount_raw_csc_indices.rds"
        )
    ),
    p = readRDS(
        file.path(
            PROJECT_DIR, "data/drop-seq", "expr_readcount_raw_csc_indptr.rds"
        )
    ),
    x = readRDS(
        file.path(
            PROJECT_DIR, "data/drop-seq", "expr_readcount_raw_csc_values.rds"
        )
    ),
    dims = readRDS(
        file.path(
            PROJECT_DIR, "data/drop-seq", "expr_readcount_raw_csc_shape.rds"
        )
    ),
    dimnames = readRDS(
        file.path(
            PROJECT_DIR, "data/drop-seq", "expr_readcount_raw_csc_dimnames.rds"
        )
    ),
    index1 = FALSE
)

dim(matrix_readcount_use)

[1] 27999 27416

rownames(matrix_readcount_use) <- paste(
    rownames(matrix_readcount_use),
    gene_symbols[rownames(matrix_readcount_use)],
    sep = "_"
)

matrix_readcount_use[1:5, 1:5] |>
    as.matrix()

                           JD126-1-2_TTTCTATATACA.bam
ENSMUSG00000051951_Xkr4                             0
ENSMUSG00000089699_Gm1992                           0
ENSMUSG00000102343_Gm37381                          0
ENSMUSG00000025900_Rp1                              0
ENSMUSG00000109048_Rp1                              0
                           JD126-1-2_CCTAGAAACCAG.bam
ENSMUSG00000051951_Xkr4                             0
ENSMUSG00000089699_Gm1992                           0
ENSMUSG00000102343_Gm37381                          0
ENSMUSG00000025900_Rp1                              0
ENSMUSG00000109048_Rp1                              0
                           JD126-1-2_TCATAGTCTATT.bam
ENSMUSG00000051951_Xkr4                             0
ENSMUSG00000089699_Gm1992                           0
ENSMUSG00000102343_Gm37381                          0
ENSMUSG00000025900_Rp1                              0
ENSMUSG00000109048_Rp1                              0
                           JD126-1-2_ATGACCTTTCCC.bam
ENSMUSG00000051951_Xkr4                             0
ENSMUSG00000089699_Gm1992                           0
ENSMUSG00000102343_Gm37381                          0
ENSMUSG00000025900_Rp1                              0
ENSMUSG00000109048_Rp1                              0
                           JD126-1-2_CATTAGTGATGG.bam
ENSMUSG00000051951_Xkr4                             0
ENSMUSG00000089699_Gm1992                           0
ENSMUSG00000102343_Gm37381                          0
ENSMUSG00000025900_Rp1                              0
ENSMUSG00000109048_Rp1                              0

Embedding

embedding <- readRDS(
    file = file.path(PROJECT_DIR, "data/drop-seq", "tsne_out_coords.rds")
)

dim(embedding)

[1] 25776     7

Check memory usage.

purrr::walk(
    list(matrix_readcount_use, embedding), \(x) {
        print(object.size(x), units = "auto", standard = "SI")
    }
)

751.2 MB
3.4 MB

Clustering

x_column <- "x"
y_column <- "y"

GEOM_POINT_SIZE <- 0.5
EMBEDDING_TITLE_PREFIX <- "t-SNE"
RASTERISED <- FALSE

Embedding

embedding |>
    tibble::rownames_to_column(var = "cell") |>
    dplyr::mutate(
        num_umis = colSums(matrix_readcount_use[, cell]),
        num_features = colSums(matrix_readcount_use[, cell] > 0),
    ) |>
    dplyr::rename(batch = batch.id) |>
    dplyr::group_by(batch) |>
    dplyr::summarise(
        num_cells = n(),
        median_umis = median(num_umis),
        median_features = median(num_features)
    ) |>
    gt::gt() |>
    gt::data_color(
        columns = c(median_umis),
        # colors = scales::col_numeric(
        #     palette = c(
        #         "green", "orange", "red"
        #     ),
        #     domain = c(500, 25000)
        colors = scales::col_numeric(
            palette = c(
                "green", "orange", "red"
            ),
            domain = NULL
        )
    ) |>
    gt::summary_rows(
        columns = c(batch),
        fns = list(
            Count = ~ n()
        ),
        decimals = 0
    ) |>
    gt::summary_rows(
        columns = c(num_cells),
        fns = list(
            Sum = ~ sum(.)
        ),
        decimals = 0
    ) |>
    gt::summary_rows(
        columns = c(num_cells:median_features),
        fns = list(
            Mean = ~ mean(.)
        ),
        decimals = 2
    ) |>
    gt::tab_header(
        title = gt::md("**Drop-Seq**; Batch")
        # subtitle = md("`gtcars` is an R dataset")
    ) |>
    gt::tab_source_note(
        source_note = gt::md("***[Here](https://www.ncbi.nlm.nih.gov/Traces/study/?query_key=5&WebEnv=MCID_61f9806c1b34bc07dc269bb6)** are the details.*")
    )

Drop-Seq; Batch
	batch	num_cells	median_umis	median_features
	JD126-1-2	335	3697.0	1839.0
	JD126A1	289	3499.0	1669.0
	JD126A5	310	3575.0	1759.5
	JD126B	625	3604.0	1733.0
	JD131-A	418	1765.0	1105.0
	JD131-B	371	1895.0	1146.0
	JD131-C	381	1894.0	1168.0
	JD135A	101	3203.0	1521.0
	JD135Y	222	1755.5	1007.0
	JD135Z	92	1890.0	966.5
	JD136A	172	2348.5	1285.5
	JD136B	269	1603.0	975.0
	JD136C	168	1140.5	734.0
	JD136DE	370	1042.5	694.0
	JD136F	251	1752.0	993.0
	JD136G	500	1341.5	827.0
	JD137A	268	2796.5	1410.5
	JD145A1	269	3332.0	1623.0
	JD145A5	267	2911.0	1503.0
	JD145B1	284	3531.5	1704.0
	JD145B5	287	1905.0	1045.0
	JD145C	169	2235.0	1056.0
	JD145D	187	1797.0	941.0
	JD145E	287	977.0	590.0
	JD145F	103	533.0	396.0
	JD146A1	530	3210.5	1555.5
	JD146A5	549	4260.0	1935.0
	JD146B1	466	4879.5	1973.0
	JD146B5	444	4035.5	1828.0
	JD146C1	452	5129.5	2004.5
	JD146C5	420	2922.0	1399.5
	JD149A	275	1118.0	699.0
	JD149B	251	823.0	531.0
	JD149C	382	1179.5	747.5
	JD149D	365	1054.0	688.0
	JD149E	236	854.5	603.0
	JD150-3A	393	1486.0	911.0
	JD150-3B	345	1624.0	1002.0
	JD150-3C	487	2110.0	1190.0
	JD150-4A	330	1652.0	962.0
	JD150-4B	370	1476.5	893.5
	JD150A	342	778.5	487.0
	JD150B	317	849.0	564.0
	JD150C	346	1372.0	888.5
	JD150D	465	1081.0	723.0
	JD164_165	272	3324.5	1935.5
	JD166_167	204	6267.5	3047.0
	JD168_169	236	4234.5	2302.5
	JD168A	252	11674.5	4225.0
	JD168B	240	4899.5	2510.5
	JD168C	182	14665.5	4678.5
	JD169A	306	7700.0	3450.0
	JD169B	178	13284.5	4515.5
	JD169C	206	10929.5	4092.0
	JD170_171	288	13620.5	4505.5
	JD170A	148	11175.0	4197.0
	JD170B	108	12311.0	4191.0
	JD170C	150	11961.5	4228.0
	JD171A	270	7281.0	3207.0
	JD171B	224	11545.0	4211.5
	JD171C	166	9291.0	3673.0
	PZ473	168	9368.5	3888.5
	PZ474	188	6470.0	3116.5
	PZ475	326	11516.5	4077.0
	PZ476	292	10177.5	3943.0
	PZ477	552	5562.5	2731.0
	PZ478	344	7171.0	3283.5
	PZ479	356	9106.0	3655.0
	PZ480	196	10180.5	3878.5
	PZ481	132	11876.5	4161.0
	PZ482	70	4286.0	2213.0
	PZ483	98	4452.5	2233.5
	PZ484	224	9902.5	3863.5
	PZ485	264	10482.5	3963.5
	PZ486	344	6220.0	2901.5
	PZ487	238	5281.5	2571.5
	PZ492	399	1618.0	960.0
	PZ493	241	1394.0	841.0
	PZ496	344	7004.5	3245.0
	PZ497	208	9547.0	3810.5
	PZ580	472	8199.5	3482.5
	PZ581	360	6468.0	3058.5
	PZ582	214	9730.5	3892.5
	PZ583	120	10060.5	3977.0
	PZ584	190	14141.5	4915.0
	PZ585	152	15586.5	4924.0
	PZ586	292	4879.0	2466.5
	PZ587	204	11968.5	4318.5
	PZ588	110	24547.5	5941.0
	PZ589	110	11439.0	3965.5
	PZ590	94	17243.5	5124.5
	PZ660	148	2893.5	1478.0
	PZ661	136	2414.5	1289.5
Count	93	—	—	—
Sum	—	25,776	—	—
Mean	—	277.16	5,734.12	2,348.55
Here are the details.

embedding <- embedding |>
    tibble::rownames_to_column(var = "cell") |>
    dplyr::mutate(
        group = dplyr::case_when(
            category %in% c("JD168") ~ "Reprogrammed",
            category %in% c("JD174") ~ "Control",
            TRUE ~ "Primary"
        ),
        group = factor(
            group,
            levels = c("Primary", "Reprogrammed", "Control")
        )
    )

p_embedding_cluster <- plot_embedding(
    data = embedding[, c(x_column, y_column)],
    color = embedding$cluster |> as.factor(),
    label = glue::glue("{EMBEDDING_TITLE_PREFIX}; Cluster"),
    color_labels = TRUE,
    color_legend = FALSE,
    sort_values = FALSE,
    rasterise = RASTERISED,
    geom_point_size = GEOM_POINT_SIZE
) +
    theme_customized_embedding()

p_embedding_UMI <- plot_embedding(
    data = embedding[, c(x_column, y_column)],
    color = log10(Matrix::colSums(matrix_readcount_use[, embedding$cell])),
    label = glue::glue("{EMBEDDING_TITLE_PREFIX}; UMI"),
    color_legend = TRUE,
    sort_values = FALSE,
    shuffle_values = TRUE,
    rasterise = RASTERISED,
    geom_point_size = GEOM_POINT_SIZE
) +
    theme_customized_embedding()

p_embedding_MT <- plot_embedding(
    data = embedding[, c(x_column, y_column)],
    color = (colSums(matrix_readcount_use[
        stringr::str_detect(
            string = stringr::str_remove(
                string = rownames(matrix_readcount_use),
                pattern = "^E.+_"
            ),
            pattern = "mt-"
        ),
    ]) / colSums(matrix_readcount_use))[embedding$cell],
    label = glue::glue("{EMBEDDING_TITLE_PREFIX}; MT %"),
    color_legend = TRUE,
    sort_values = TRUE,
    shuffle_values = FALSE,
    rasterise = RASTERISED,
    geom_point_size = GEOM_POINT_SIZE
) + theme_customized_embedding()

p_embedding_group <- plot_embedding(
    data = embedding[, c(x_column, y_column)],
    color = embedding$group |> as.factor(),
    label = glue::glue("{EMBEDDING_TITLE_PREFIX}; Group"),
    color_labels = FALSE,
    color_legend = TRUE,
    sort_values = FALSE,
    shuffle_values = TRUE,
    rasterise = RASTERISED,
    geom_point_size = GEOM_POINT_SIZE
) +
    theme_customized_embedding() +
    ggplot2::scale_color_manual(
        values = c(
            Primary = "#00AFBB",
            Reprogrammed = "#8BC34A",
            Control = "#E7B800"
        )
    )

embedding |>
    dplyr::mutate(
        num_umis = colSums(matrix_readcount_use[, cell]),
        num_features = colSums(matrix_readcount_use[, cell] > 0),
    ) |>
    dplyr::group_by(cluster) |>
    dplyr::summarise(
        num_cells = n(),
        median_umis = median(num_umis),
        median_features = median(num_features)
    ) |>
    gt::gt() |>
    gt::data_color(
        columns = c(median_umis),
        colors = scales::col_numeric(
            palette = c(
                "green", "orange", "red"
            ),
            domain = NULL
        )
    ) |>
    gt::summary_rows(
        columns = c(cluster),
        fns = list(
            Count = ~ n()
        ),
        decimals = 0
    ) |>
    gt::summary_rows(
        columns = c(num_cells),
        fns = list(
            Sum = ~ sum(.)
        ),
        decimals = 0
    ) |>
    gt::tab_header(
        title = gt::md("**Drop-Seq**; Clustering")
    )

Drop-Seq; Clustering
	cluster	num_cells	median_umis	median_features
	1	2875	13196.0	4578.0
	2	2798	8594.0	3639.5
	3	2786	6165.0	2905.0
	4	2710	933.0	561.5
	5	2356	1505.5	1013.0
	6	2139	2590.0	1249.0
	7	1934	2238.5	1347.0
	8	1575	5225.0	1891.0
	9	1075	5719.0	2203.0
	10	871	12955.0	4294.0
	11	708	4111.0	1872.5
	12	681	3495.0	1956.0
	13	669	2438.0	1483.0
	14	556	725.0	530.5
	15	478	2661.0	1606.5
	16	474	1368.0	906.5
	17	254	4995.5	2172.5
	18	236	1737.0	256.5
	19	181	620.0	372.0
	20	165	9490.0	3679.0
	21	133	5234.0	2464.0
	22	122	6274.0	2877.0
Count	22	—	—	—
Sum	—	25,776	—	—

purrr::reduce(
    list(
        p_embedding_cluster,
        p_embedding_UMI,
        p_embedding_MT,
        p_embedding_group
    ),
    `+`
) +
    patchwork::plot_layout(ncol = 2) +
    patchwork::plot_annotation(
        theme = ggplot2::theme(plot.margin = ggplot2::margin())
    )

library(formattable)

Attaching package: 'formattable'

The following object is masked from 'package:patchwork':

    area

embedding |>
    dplyr::mutate(
        num_umis = colSums(matrix_readcount_use[, cell]),
        num_features = colSums(matrix_readcount_use[, cell] > 0),
    ) |>
    dplyr::group_by(group) |>
    dplyr::summarise(
        num_cells = n(),
        median_umis = median(num_umis),
        median_features = median(num_features)
    ) |>
    formattable::formattable(
        list(
            # num_cells = formattable::color_tile("transparent", "lightpink"),
            num_cells = formattable::color_bar("Lightpink"),
            median_umis = formattable::color_bar("lightgreen"),
            median_features = formattable::color_bar("lightblue")
        ),
        full_width = FALSE,
        caption = "Drop-Seq; Group"
    )

Drop-Seq; Group

group	num_cells	median_umis	median_features
Primary	15684	2030	1139
Reprogrammed	8730	8532	3562
Control	1362	9383	3789

purrr::map(levels(embedding$group), \(x) {
    plot_embedding(
        data = embedding[, c(x_column, y_column)],
        color = as.integer(embedding$group == x) |> as.factor(),
        label = glue::glue(
            "{EMBEDDING_TITLE_PREFIX}; {x}: {sum(embedding$group == x)}"
        ),
        color_labels = FALSE,
        color_legend = FALSE,
        sort_values = TRUE,
        shuffle_values = FALSE,
        rasterise = RASTERISED,
        geom_point_size = GEOM_POINT_SIZE
    ) +
        theme_customized_embedding() +
        ggplot2::scale_color_manual(
            values = c("grey70", "salmon")
        )
}) |>
    purrr::reduce(`+`) +
    patchwork::plot_layout(ncol = 3) +
    patchwork::plot_annotation(
        theme = ggplot2::theme(plot.margin = ggplot2::margin())
    )

Extract colors from the initial plots to keep colors consistent.

color_palette <- ggplot2::ggplot_build(p_embedding_cluster)$data[[1]] |>
    dplyr::select(color = colour, cluster = group) |>
    unique() |>
    dplyr::arrange(cluster)

color_palette <- setNames(
    object = color_palette$color,
    nm = color_palette$cluster
)

scales::show_col(color_palette, borders = NA)

cell_type_segments <- data.frame(
    x = c(-40, -40, -40, -40, -69, -35, 40, 40, 40, 65, 28, -12, -57.5, -57.5, -57.5),
    y = c(60, 60, 60, 60, 20, -55, -70, -70, -70, 50, 55, 65, -52.5, -52.5, -52.5),
    xend = c(-30, -44, -60, -30, -69, -35, 32, 0, 40, 65, 35, -20, -17.5, -53, -37.5),
    yend = c(35, 15, -5, -22, 35, -65, -58, -66, -45, 36, 60, 70, -25, -45, -47.5),
    cluster = c(4, 6, 9, 11, 8, 12, 7, 13, 14, 16, 5, 15, 2, 20, 12)
)

cell_type_labels <- data.frame(
    x = c(-40, -61.5, -43, 52, 56, 46, -22, -58),
    y = c(65, 45, -70, -70, 55, 64, 75, -58),
    label = c(
        "Atrial myocytes",
        "Ventricular\nmyocytes",
        "Epicardial cells",
        "Fibroblasts",
        "Lymphocytes",
        "Endothelial cells",
        "Macrophages",
        "MEF-derived"
    )
)

cluster_labels <- embedding |>
    dplyr::group_by(cluster) |>
    dplyr::summarise(
        x = median(x),
        y = median(y)
    ) |>
    as.data.frame()

plot_embedding(
    data = embedding[, c(x_column, y_column)],
    color = embedding |>
        dplyr::mutate(
            color_group = dplyr::case_when(
                group %in% c("Reprogrammed", "Control") ~ "MEF-derived",
                TRUE ~ as.character(cluster)
            )
        ) |>
        dplyr::pull(color_group) |> as.factor(),
    label = NULL,
    color_labels = FALSE,
    color_legend = FALSE,
    sort_values = FALSE,
    rasterise = RASTERISED,
    geom_point_size = GEOM_POINT_SIZE
) +
    theme_customized_embedding(void = TRUE) +
    ggplot2::scale_color_manual(
        values = c(color_palette, "MEF-derived" = "salmon")
    ) +
    ggplot2::geom_segment(
        data = cell_type_segments,
        ggplot2::aes(x = x, xend = xend, y = y, yend = yend),
        color = "grey50",
        size = .2
    ) +
    ggplot2::geom_text(
        data = cell_type_labels,
        ggplot2::aes(x, y, label = label),
        color = c(rep("black", 2), "#00BFC4", rep("black", 4), "#FF5722"),
        size = 2.8,
        family = "Arial"
    ) +
    ggplot2::annotate(
        geom = "text",
        family = "Arial",
        x = cluster_labels[, "x"],
        y = cluster_labels[, "y"], label = cluster_labels[, 1],
        parse = TRUE,
        size = 2,
        color = c("black")
    )

pdf_width <- 104
pdf_height <- 74
file_name <- glue::glue(
    "Rplot_embedding_dropseq_{EMBEDDING_TITLE_PREFIX}_",
    "cell_group_{pdf_width}_{pdf_height}.pdf"
)
if (!file.exists(file_name)) {
    ggplot2::ggsave(
        filename = file_name,
        useDingbats = FALSE,
        plot = ggplot2::last_plot(),
        device = NULL,
        path = NULL,
        scale = 1,
        width = pdf_width,
        height = pdf_height,
        units = c("mm"),
    )
}

Composition

Bar charts indicating the cellar composition of t-SNE clusters defined in Figure 1B.

calc_group_composition(
    data = embedding,
    x = "cluster",
    group = "group"
) |>
    dplyr::mutate(
        cluster = factor(
            cluster,
            levels = c(
                19, 18,
                9, 6, 11, 4, 8, 21,
                3, 1, 2, 10, 17, 22, 20,
                13, 7, 14, 12, 15, 5, 16
            )
        )
    ) |>
    plot_barplot(
        x = "cluster",
        y = "percentage",
        z = "group",
        legend_ncol = 1
    ) +
    ggplot2::scale_fill_manual(
        values = c(
            Primary = "#00AFBB",
            Reprogrammed = "#8BC34A",
            Control = "#E7B800"
        )
    )

Expression

Embedding

FEATURES_SELECTED <- c(
    "ENSMUSG00000009471_Myod1",
    "ENSMUSG00000026414_Tnnt2",
    "ENSMUSG00000016458_Wt1",
    "ENSMUSG00000025105_Bnc1",
    "ENSMUSG00000049382_Krt8",
    "ENSMUSG00000079018_Ly6c1",
    "ENSMUSG00000049436_Upk1b",
    "ENSMUSG00000021391_Cenpp"
)

purrr::map(FEATURES_SELECTED, \(x) {
    selected_feature <- x

    cat(selected_feature, "\n")
    values <- log10(
        calc_cpm(matrix_readcount_use[, embedding$cell])
        [selected_feature, ] + 1
    )

    p1 <- plot_embedding(
        data = embedding[, c(x_column, y_column)],
        color = values,
        label = paste(
            EMBEDDING_TITLE_PREFIX,
            selected_feature |> stringr::str_remove(pattern = "^E.+_"),
            sep = "; "
        ),
        color_legend = TRUE,
        sort_values = TRUE,
        rasterise = RASTERISED,
        geom_point_size = GEOM_POINT_SIZE * 1.25,
        na_value = "grey80"
    ) +
        theme_customized_embedding()

    return(p1)
}) |>
    # unlist(recursive = FALSE) |>
    purrr::reduce(`+`) +
    patchwork::plot_layout(ncol = 2, byrow = FALSE) +
    patchwork::plot_annotation(
        theme = ggplot2::theme(plot.margin = ggplot2::margin())
    )

ENSMUSG00000009471_Myod1 
ENSMUSG00000026414_Tnnt2 
ENSMUSG00000016458_Wt1 
ENSMUSG00000025105_Bnc1 
ENSMUSG00000049382_Krt8 
ENSMUSG00000079018_Ly6c1 
ENSMUSG00000049436_Upk1b 
ENSMUSG00000021391_Cenpp 

Violin

Violin plot illustrating the expression of cardiac markers from single-cell expression data derived from P0 mouse heart and reprogrammed/uninfected MEFs.

labels_y <- c(
    19, 18,
    9, 6, 11, 4, 8, 21,
    3, 1, 2, 10, 17, 22, 20,
    13, 7, 14, 12, 15, 5, 16
)

plot_violin(
    cells = embedding |>
        dplyr::mutate(
            cluster = factor(
                cluster,
                levels = labels_y
            )
        ) |>
        split(~cluster) |>
        purrr::map(\(x) {
            x |> dplyr::pull(cell)
        }),
    features = c(
        "ENSMUSG00000026414_Tnnt2",
        "ENSMUSG00000013936_Myl2",
        "ENSMUSG00000042045_Sln",
        "ENSMUSG00000001506_Col1a1",
        "ENSMUSG00000005836_Gata6",
        "ENSMUSG00000005583_Mef2c"
    ),
    matrix_cpm = calc_cpm(matrix_readcount_use)
) +
    theme_customized_violin(
        axis_text_color_y = rev(color_palette[as.character(labels_y)])
    ) +
    ggplot2::scale_fill_manual(
        values = color_palette[as.character(labels_y)]
    ) +
    ggplot2::scale_color_manual(
        values = color_palette[as.character(labels_y)]
    )

Bar plot

barplot_helper <- function(cells, features, matrix_readcount) {
    purrr::map(names(cells), \(x) {
        calc_cpm(matrix_readcount)[
            features,
            colnames(matrix_readcount) %in% cells[[x]]
        ] |>
            Matrix::rowMeans() |>
            tibble::enframe(name = "feature") |>
            dplyr::mutate(group = x)
    }) |>
        dplyr::bind_rows()
}

Expression of Myod1 in MEF-derived cells.

cells_barplot <- list(
    embedding |>
        dplyr::filter(
            cluster == 20 & category %in% ("JD168")
        ) |>
        dplyr::pull(cell),
    embedding |>
        dplyr::filter(
            cluster != 20 & category %in% ("JD168")
        ) |>
        dplyr::pull(cell),
    embedding |>
        dplyr::filter(
            category %in% ("JD174")
        ) |>
        dplyr::pull(cell)
)
names(cells_barplot) <- c("Cluster 20", "Other", "Uninfected")

features_barplot <- c(
    "ENSMUSG00000009471_Myod1",
    "ENSMUSG00000016458_Wt1",
    "ENSMUSG00000025105_Bnc1"
)

barplot_helper(cells_barplot, features_barplot, matrix_readcount_use) |>
    dplyr::filter(feature == "ENSMUSG00000009471_Myod1") |>
    dplyr::mutate(
        feature = stringr::str_remove(string = feature, pattern = "^.+_")
    ) |>
    dplyr::mutate(value = log10(value + 1)) |>
    plot_barplot_simple(
        x = "group",
        y = "value",
        z = "feature",
        y_title = expression("Avg expr; log"[10] * " (CPM + 1)")
    ) +
    theme_customized_violin(
        strip_background_fill = "grey80",
        panel_border_color = "black",
        axis_text_x_angle = c(90, 1, 0.5)
    ) +
    ggplot2::scale_fill_manual(
        values = c(
            c(
                "#FF5722",
                "grey35",
                "grey35"
            )
        )
    )

Expression of Wt1 and Bnc1 in MEF-derived cells.

cells_barplot <- list(
    embedding |>
        dplyr::filter(
            cluster == 12,
            !category %in% c("JD168", "JD174")
        ) |>
        dplyr::pull(cell),
    embedding |>
        dplyr::filter(
            cluster == 12 & category %in% ("JD168")
        ) |>
        dplyr::pull(cell),
    embedding |>
        dplyr::filter(
            cluster != 12 & category %in% ("JD168")
        ) |>
        dplyr::pull(cell),
    embedding |>
        dplyr::filter(
            category %in% ("JD174")
        ) |>
        dplyr::pull(cell)
)
names(cells_barplot) <- c("Primarry epi.", "Cluster 12", "Other", "Uninfected")

barplot_helper(cells_barplot, features_barplot, matrix_readcount_use) |>
    dplyr::filter(feature != "ENSMUSG00000009471_Myod1") |>
    dplyr::mutate(
        feature = stringr::str_remove(string = feature, pattern = "^.+_"),
        feature = factor(
            feature,
            levels = c("Wt1", "Bnc1")
        )
    ) |>
    dplyr::mutate(
        value = log10(value + 1),
        group = factor(
            group,
            levels = names(cells_barplot)
        )
    ) |>
    plot_barplot_simple(
        x = "group",
        y = "value",
        z = "feature",
        y_title = expression("Avg expr; log"[10] * " (CPM + 1)")
    ) +
    theme_customized_violin(
        strip_background_fill = "grey80",
        panel_border_color = "black",
        axis_text_x_angle = c(90, 1, 0.5)
    ) +
    ggplot2::scale_fill_manual(
        values = c(
            c(
                "#00BFC4",
                "#FF5722",
                "grey35",
                "grey35"
            )
        )
    )

Lollipop

10 TFs differentially expressed in primary epicardial cells compared with uninfected MEFs. Shown is the expression of 10F in MEFs and P0 mouse heart cells.

cells_lollipop <- embedding |>
    dplyr::filter(
        !cluster %in% c(1, 2, 3, 10, 17, 20, 22),
        !category %in% c("JD168", "JD174")
    ) |>
    split(~cluster) |>
    purrr::map(\(x) {
        x |> dplyr::pull(cell)
    })

cells_lollipop <- cells_lollipop[
    purrr::map_lgl(cells_lollipop, \(x) length(x) > 0)
]

cells_lollipop$control <- embedding |>
    dplyr::filter(category %in% c("JD168", "JD174")) |>
    dplyr::pull(cell)


labels_y <- c(
    "12", "control", "19", "18", "9", "6", "11", "4", "8",
    "21", "13", "7", "14", "15", "5", "16"
)

plot_lollipop(
    cells = cells_lollipop[labels_y],
    features = c(
        "ENSMUSG00000026628_Atf3",
        "ENSMUSG00000016458_Wt1",
        "ENSMUSG00000025105_Bnc1",
        "ENSMUSG00000051910_Sox6",
        "ENSMUSG00000045680_Tcf21",
        "ENSMUSG00000038193_Hand2",
        "ENSMUSG00000031965_Tbx20",
        "ENSMUSG00000032419_Tbx18",
        "ENSMUSG00000005836_Gata6",
        "ENSMUSG00000036098_Myrf"
    ),
    matrix_cpm = calc_cpm(matrix_readcount_use)
) +
    ggplot2::scale_y_discrete(
        name = NULL,
        labels = c(
            "Epicardial cell",
            "Uninfected MEF",
            "CM 2",
            "CM 1",
            "Atrial CM",
            "Atrial CM",
            "Atrial CM",
            "Atrial CM",
            "Ventricular CM",
            "CM 3",
            "Cardiac fibroblast",
            "Cardiac fibroblast",
            "Cardiac fibroblast",
            "Macrophage",
            "Endothelial cell",
            "Lymphocyte"
        ) |> rev()
    ) +
    theme_customized_violin(
        axis_text_x_angle = c(-45, 1, 0.5),
        axis_text_color_y = c(color_palette, control = "grey35")[rev(labels_y)],
        panel_grid_major = TRUE
    ) %+replace%
    ggplot2::theme(
        legend.margin = ggplot2::margin(
            t = 0, r = 0, b = 0, l = 0, unit = "mm"
        ),
        legend.key.size = ggplot2::unit(2.5, "mm"),
        legend.key.width = ggplot2::unit(4.0, "mm"),
        legend.text = ggplot2::element_text(family = "Arial", size = 5),
        legend.title = ggplot2::element_text(family = "Arial", size = 6),
        legend.position = "bottom",
        legend.box = "vertical" # "horizontal"
    ) +
    ggplot2::scale_size(
        name = "% cells expressing gene",
        breaks = seq(0, 1, .2),
        labels = seq(0, 1, .2) * 100,
        limits = c(0, 1),
        range = c(0, 6),
        guide = ggplot2::guide_legend(
            title.position = "top",
            title.hjust = 0.5,
            label.position = "bottom",
            nrow = 1,
            byrow = TRUE,
            order = 1
        )
    ) +
    ggplot2::scale_color_viridis_c(
        name = expression(paste("Avg expr; log"[10], "( CPM + 1)")),
        guide = ggplot2::guide_colourbar(
            title.position = "top",
            title.hjust = 1,
            barwidth = 5,
            barheight = 0.6,
            direction = "horizontal",
            order = 2
        )
    )

Gene Ontology enrichment

Gene ontology analysis for genes highly expressed in Cluster 20.

suppressPackageStartupMessages(library(topGO))

groupGOTerms:   GOBPTerm, GOMFTerm, GOCCTerm environments built.

packageVersion("topGO")

[1] '2.46.0'

packageVersion("org.Mm.eg.db")

[1] '3.14.0'

de_paired <- detect_de(
    cell_group_a = embedding |>
        dplyr::filter(
            category == "JD168",
            cluster == 20
        ) |>
        dplyr::pull(cell),
    cell_group_b = embedding |>
        dplyr::filter(
            category == "JD168",
            cluster != 20
        ) |>
        dplyr::pull(cell),
    matrix_readcount = matrix_readcount_use,
    matrix_cpm = calc_cpm(matrix_readcount_use)
    # only_enrichment = TRUE
)

de_paired |> head()

                           log2_effect          pval positive_frac_a
ENSMUSG00000037139_Myom3      8.166286 3.152517e-219           0.604
ENSMUSG00000079588_Tmem182    8.050808 5.009801e-140           0.402
ENSMUSG00000026251_Chrnd      7.992256  5.311037e-81           0.238
ENSMUSG00000087591_Gm14635    7.975343 2.164927e-110           0.323
ENSMUSG00000101680_Gm29015    7.963194  4.971533e-89           0.262
ENSMUSG00000102717_Gm37759    7.876779 1.253173e-149           0.439
                           positive_frac_b norm_reads_mean_a norm_reads_mean_b
ENSMUSG00000037139_Myom3             0.002         1.5196756      0.0018586658
ENSMUSG00000079588_Tmem182           0.002         0.8570739      0.0008143528
ENSMUSG00000026251_Chrnd             0.001         0.2803849      0.0004863366
ENSMUSG00000087591_Gm14635           0.001         0.6266815      0.0006813209
ENSMUSG00000101680_Gm29015           0.001         0.4385369      0.0020254669
ENSMUSG00000102717_Gm37759           0.002         0.7279072      0.0018620079
                           log2_fc_norm_reads cpm_meam_a cpm_meam_b log2_fc_cpm
ENSMUSG00000037139_Myom3             7.011140  177.63879 0.21767902    9.607811
ENSMUSG00000079588_Tmem182           6.325136  100.13546 0.09540360    9.891957
ENSMUSG00000026251_Chrnd             4.791384   32.73270 0.05698473    8.933125
ENSMUSG00000087591_Gm14635           5.897410   73.21501 0.07976347    9.671992
ENSMUSG00000101680_Gm29015           5.221062   51.29588 0.23699159    7.698518
ENSMUSG00000102717_Gm37759           5.959019   85.12481 0.21819530    8.543336
                                pval_adj
ENSMUSG00000037139_Myom3   6.067921e-219
ENSMUSG00000079588_Tmem182 8.772996e-140
ENSMUSG00000026251_Chrnd    8.352273e-81
ENSMUSG00000087591_Gm14635 3.606077e-110
ENSMUSG00000101680_Gm29015  7.951628e-89
ENSMUSG00000102717_Gm37759 2.229182e-149

dim(de_paired)

[1] 1126   11

genes_of_interest <- rownames(subset(de_paired, log2_effect > 0))
gene_universe <- rownames(matrix_readcount_use)
genes_formatted <- factor(as.integer(gene_universe %in% genes_of_interest))
names(genes_formatted) <- gene_universe
names(genes_formatted) <- names(genes_formatted) |>
    stringr::str_remove(pattern = "_.+$")

topgo_data <- new(
    "topGOdata",
    ontology = "BP",
    allGenes = genes_formatted,
    annot = annFUN.org,
    mapping = "org.Mm.eg.db",
    ID = "Ensembl"
)
## 
## Building most specific GOs .....
## Loading required package: org.Mm.eg.db
## 
##  ( 12529 GO terms found. )
## 
## Build GO DAG topology ..........
##  ( 15987 GO terms and 36498 relations. )
## 
## Annotating nodes ...............
##  ( 21119 genes annotated to the GO terms. )

topgo_out_classic_fisher <- topGO::runTest(
    topgo_data,
    algorithm = "classic",
    statistic = "fisher"
)
## 
##           -- Classic Algorithm -- 
## 
##       the algorithm is scoring 5197 nontrivial nodes
##       parameters: 
##           test statistic: fisher

NUM_GO_TERMS <- 15

# prepare data
enriched_gos <- topGO::GenTable(topgo_data,
    classicFisher = topgo_out_classic_fisher,
    topNodes = NUM_GO_TERMS
) |>
    dplyr::mutate(
        classicFisher = as.numeric(
            stringr::str_replace(classicFisher, "< ", "")
        ),
        Term = factor(Term(GO.ID),
            levels = rev(Term(GO.ID))
        )
    )

ggplot2::ggplot(
    enriched_gos,
    ggplot2::aes(
        y = -log10(classicFisher),
        x = Term
    )
) +
    ggplot2::geom_bar(stat = "identity", fill = "grey70") +
    ggplot2::coord_flip() +
    ggplot2::labs(x = NULL) +
    ggplot2::theme_classic() +
    ggplot2::scale_y_continuous(
        name = expression(paste("-log"[10], " (p-value)")),
        breaks = c(0, 30),
        labels = scales::math_format(10^.x)
    ) +
    ggplot2::annotate(
        geom = "text",
        x = seq_len(nrow(enriched_gos)),
        y = rep(1, nrow(enriched_gos)),
        label = rev(enriched_gos$Term),
        # vjust = "inward",
        hjust = "inward",
        size = 2.5,
        family = "Arial"
    ) +
    ggplot2::theme(
        axis.title = ggplot2::element_text(family = "Arial", size = 8),
        axis.title.y = ggplot2::element_blank(),
        axis.text = ggplot2::element_text(family = "Arial", size = 7),
        axis.text.y = ggplot2::element_blank(),
        axis.ticks.y = ggplot2::element_blank(),
        axis.line.y = ggplot2::element_blank(),
        legend.text = ggplot2::element_text(family = "Arial", size = 8),
        legend.title = ggplot2::element_text(family = "Arial", size = 8)
    )

Enrichment of exogenous factors

# clusters_selected <- c(1, 2, 3, 10, 12, 17, 20, 22)

features_selected_43 <- c(
    "ENSMUSG00000025930_Msc",
    "ENSMUSG00000026313_Hdac4",
    "ENSMUSG00000026565_Pou2f1",
    "ENSMUSG00000026923_Notch1",
    "ENSMUSG00000015846_Rxra",
    "ENSMUSG00000015627_Gata5",
    "ENSMUSG00000025860_Xiap",
    "ENSMUSG00000040289_Hey1",
    "ENSMUSG00000027833_Shox2",
    "ENSMUSG00000001419_Mef2d",
    "ENSMUSG00000028800_Hdac1",
    "ENSMUSG00000086369_E330017L17Rik",
    "ENSMUSG00000028949_Smarcd3",
    "ENSMUSG00000048450_Msx1",
    "ENSMUSG00000042002_Foxn4",
    "ENSMUSG00000018604_Tbx3",
    "ENSMUSG00000018263_Tbx5",
    "ENSMUSG00000063568_Jazf1",
    "ENSMUSG00000009471_Myod1",
    "ENSMUSG00000030557_Mef2a",
    "ENSMUSG00000030551_Nr2f2",
    "ENSMUSG00000030544_Mesp1",
    "ENSMUSG00000019789_Hey2",
    "ENSMUSG00000019777_Hdac2",
    "ENSMUSG00000020167_Tcf3",
    "ENSMUSG00000038193_Hand2",
    "ENSMUSG00000079033_Mef2b",
    "ENSMUSG00000021944_Gata4",
    "ENSMUSG00000032419_Tbx18",
    "ENSMUSG00000020160_Meis1",
    "ENSMUSG00000037335_Hand1",
    "ENSMUSG00000020542_Myocd",
    "ENSMUSG00000000093_Tbx2",
    "ENSMUSG00000021469_Msx2",
    "ENSMUSG00000005583_Mef2c",
    "ENSMUSG00000042258_Isl1",
    "ENSMUSG00000009739_Pou6f1",
    "ENSMUSG00000001288_Rarg",
    "ENSMUSG00000015579_Nkx2-5",
    "ENSMUSG00000023067_Cdkn1a",
    "ENSMUSG00000024063_Lbh",
    "ENSMUSG00000005836_Gata6",
    "ENSMUSG00000024515_Smad4"
)

clusters_selected <- c(20, 12, 2)

enriched_factors <- do.call(
    rbind.data.frame,
    lapply(clusters_selected, \(x) {
        cells_1 <- embedding$cell[
            embedding$cluster == x & embedding$category == "JD168"
        ]
        cells_2 <- embedding$cell[
            embedding$cluster != x & embedding$category == "JD168"
        ]
        cat(x, length(cells_1), length(cells_2), "\n")

        de_paired <- detect_de(
            cell_group_a = cells_1,
            cell_group_b = cells_2,
            matrix_readcount = matrix_readcount_use,
            matrix_cpm = calc_cpm(matrix_readcount_use),
            only_enrichment = TRUE
        ) |>
            dplyr::mutate(category = x) |>
            tibble::rownames_to_column(var = "feature") |>
            dplyr::filter(feature %in% features_selected_43)
    })
) |>
    dplyr::filter(category %in% clusters_selected) |>
    dplyr::mutate(
        category = factor(category,
            levels = clusters_selected
        ),
        symbol = stringr::str_remove(
            string = feature,
            pattern = "^.+_"
        )
    )

20 164 8566 
12 201 8529 
2 2445 6285 

Differential expression analysis of 48F in MEF-derived clusters, as compared with all other reprogrammed cells. Each dot represents a gene (colored by fold change and sized by p value).

`%+replace%` <- ggplot2::`%+replace%`

ggplot2::ggplot() +
    ggplot2::geom_abline(intercept = 0, slope = 1, linetype = 2) +
    ggplot2::geom_point(
        data = enriched_factors,
        ggplot2::aes(positive_frac_b,
            positive_frac_a,
            size = -log10(pval_adj),
            color = log2_effect
        ),
        alpha = .8,
        stroke = 0, shape = 16
    ) +
    ggplot2::facet_wrap(
        ~category,
        nrow = 1,
        labeller = ggplot2::labeller(
            category = setNames(
                object = paste("Cluster", clusters_selected),
                nm = clusters_selected
            )
        )
    ) +
    ggplot2::coord_fixed() +
    ggplot2::scale_color_viridis_c(
        name = expression(paste("Log"[2], " effect"))
    ) +
    ggplot2::scale_size_continuous(
        name = expression(paste("-log"[10], " (p-value)"))
    ) +
    ggplot2::guides(
        color = ggplot2::guide_colorbar(order = 1),
        size = ggplot2::guide_legend(order = 2)
    ) +
    ggplot2::scale_x_continuous(
        name = "Expr (%, other reprogrammed cells)",
        limits = c(0, 1), breaks = seq(0, 1, .2)
    ) +
    ggplot2::scale_y_continuous(
        name = "Expr (%, indicated cluster)",
        limits = c(0, 1), breaks = seq(0, 1, .2)
    ) +
    theme_customized_violin() %+replace%
    ggplot2::theme(
        legend.background = ggplot2::element_blank(),
        legend.margin = ggplot2::margin(t = 0, r = 0, b = 0, l = 0, unit = "mm"),
        legend.key.size = ggplot2::unit(2.5, "mm"),
        legend.text = ggplot2::element_text(family = "Arial", size = 6),
        legend.title = ggplot2::element_text(family = "Arial", size = 7),
        legend.position = "bottom",
        legend.box = "horizontal",
        legend.box.background = ggplot2::element_blank()
    ) +
    ggrepel::geom_text_repel(
        data = subset(enriched_factors, symbol %in% c(
            "Myod1", "Mef2b",
            "Mef2c", "Hand2", "Gata6"
        )),
        ggplot2::aes(
            positive_frac_b,
            positive_frac_a,
            label = symbol
        ),
        #
        size = 5 / ggplot2::.pt,
        family = "Arial",
        box.padding = .2,
        point.padding = .2,
        nudge_y = .15,
        arrow = ggplot2::arrow(length = ggplot2::unit(.02, "npc")),
        segment.color = "grey35",
        color = "black"
    )

R session info

devtools::session_info()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.1.2 (2021-11-01)
 os       macOS Monterey 12.2
 system   aarch64, darwin20.6.0
 ui       unknown
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Chicago
 date     2022-02-01
 pandoc   2.17.1.1 @ /opt/homebrew/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package          * version     date (UTC) lib source
 AnnotationDbi    * 1.56.2      2021-11-09 [1] Bioconductor
 assertthat         0.2.1       2019-03-21 [1] CRAN (R 4.1.1)
 backports          1.4.1       2021-12-13 [1] CRAN (R 4.1.2)
 Biobase          * 2.54.0      2021-10-26 [1] Bioconductor
 BiocGenerics     * 0.40.0      2021-10-26 [1] Bioconductor
 Biostrings         2.62.0      2021-10-26 [1] Bioconductor
 bit                4.0.4       2020-08-04 [1] CRAN (R 4.1.1)
 bit64              4.0.5       2020-08-30 [1] CRAN (R 4.1.1)
 bitops             1.0-7       2021-04-24 [1] CRAN (R 4.1.1)
 blob               1.2.2       2021-07-23 [1] CRAN (R 4.1.1)
 brio               1.1.3       2021-11-30 [1] CRAN (R 4.1.2)
 broom              0.7.12      2022-01-28 [1] CRAN (R 4.1.2)
 cachem             1.0.6       2021-08-19 [1] CRAN (R 4.1.1)
 callr              3.7.0       2021-04-20 [1] CRAN (R 4.1.1)
 cellranger         1.1.0       2016-07-27 [1] CRAN (R 4.1.1)
 checkmate          2.0.0       2020-02-06 [1] CRAN (R 4.1.1)
 cli                3.1.1       2022-01-20 [1] CRAN (R 4.1.2)
 colorspace         2.0-2       2021-06-24 [1] CRAN (R 4.1.1)
 commonmark         1.7         2018-12-01 [1] CRAN (R 4.1.1)
 crayon             1.4.2       2021-10-29 [1] CRAN (R 4.1.1)
 data.table         1.14.2      2021-09-27 [1] CRAN (R 4.1.1)
 DBI                1.1.2       2021-12-20 [1] CRAN (R 4.1.2)
 dbplyr             2.1.1       2021-04-06 [1] CRAN (R 4.1.1)
 desc               1.4.0       2021-09-28 [1] CRAN (R 4.1.1)
 devtools           2.4.3.9000  2022-01-27 [1] Github (r-lib/devtools@41280ac)
 digest             0.6.29      2021-12-01 [1] CRAN (R 4.1.2)
 dplyr            * 1.0.7.9000  2022-01-30 [1] Github (tidyverse/dplyr@0d8fb82)
 dtplyr             1.2.1.9000  2022-01-27 [1] Github (tidyverse/dtplyr@8542cf1)
 ellipsis           0.3.2       2021-04-29 [1] CRAN (R 4.1.1)
 evaluate           0.14        2019-05-28 [1] CRAN (R 4.1.1)
 extrafont        * 0.17        2014-12-08 [1] CRAN (R 4.1.1)
 extrafontdb        1.0         2012-06-11 [1] CRAN (R 4.1.1)
 fansi              1.0.2       2022-01-14 [1] CRAN (R 4.1.2)
 farver             2.1.0       2021-02-28 [1] CRAN (R 4.1.1)
 fastmap            1.1.0       2021-01-25 [1] CRAN (R 4.1.1)
 forcats          * 0.5.1.9000  2022-01-27 [1] Github (tidyverse/forcats@b4dade0)
 formattable      * 0.2.1       2021-01-07 [1] CRAN (R 4.1.2)
 fs                 1.5.2.9000  2022-01-27 [1] Github (r-lib/fs@6d1182f)
 gargle             1.2.0       2021-07-02 [1] CRAN (R 4.1.1)
 generics           0.1.1       2021-10-25 [1] CRAN (R 4.1.1)
 GenomeInfoDb       1.30.0      2021-10-26 [1] Bioconductor
 GenomeInfoDbData   1.2.7       2022-01-06 [1] Bioconductor
 ggplot2          * 3.3.5.9000  2022-01-27 [1] Github (tidyverse/ggplot2@c89c265)
 ggrepel            0.9.1       2021-01-15 [1] CRAN (R 4.1.1)
 glue               1.6.1.9000  2022-01-27 [1] Github (tidyverse/glue@3da70df)
 GO.db            * 3.14.0      2022-01-27 [1] Bioconductor
 googledrive        2.0.0       2021-07-08 [1] CRAN (R 4.1.1)
 googlesheets4      1.0.0       2021-07-21 [1] CRAN (R 4.1.1)
 graph            * 1.72.0      2021-10-26 [1] Bioconductor
 gt                 0.3.1       2021-08-07 [1] CRAN (R 4.1.2)
 gtable             0.3.0.9000  2022-01-27 [1] Github (r-lib/gtable@a0bd272)
 haven              2.4.3       2021-08-04 [1] CRAN (R 4.1.1)
 hms                1.1.1       2021-09-26 [1] CRAN (R 4.1.1)
 htmltools          0.5.2       2021-08-25 [1] CRAN (R 4.1.1)
 htmlwidgets        1.5.4       2021-09-08 [1] CRAN (R 4.1.1)
 httr               1.4.2       2020-07-20 [1] CRAN (R 4.1.1)
 IRanges          * 2.28.0      2021-10-26 [1] Bioconductor
 jsonlite           1.7.3       2022-01-17 [1] CRAN (R 4.1.2)
 KEGGREST           1.34.0      2021-10-26 [1] Bioconductor
 knitr              1.37        2021-12-16 [1] CRAN (R 4.1.2)
 labeling           0.4.2       2020-10-20 [1] CRAN (R 4.1.1)
 lattice            0.20-45     2021-09-22 [2] CRAN (R 4.1.2)
 lifecycle          1.0.1       2021-09-24 [1] CRAN (R 4.1.1)
 lubridate          1.8.0       2022-01-27 [1] Github (tidyverse/lubridate@566590f)
 magrittr           2.0.2       2022-01-26 [1] CRAN (R 4.1.2)
 Matrix           * 1.4-0       2021-12-08 [2] CRAN (R 4.1.2)
 matrixStats        0.61.0      2021-09-17 [1] CRAN (R 4.1.1)
 memoise            2.0.1       2021-11-26 [1] CRAN (R 4.1.2)
 modelr             0.1.8.9000  2022-01-27 [1] Github (tidyverse/modelr@16168e0)
 munsell            0.5.0       2018-06-12 [1] CRAN (R 4.1.1)
 org.Mm.eg.db     * 3.14.0      2022-01-27 [1] Bioconductor
 patchwork        * 1.1.0.9000  2022-01-27 [1] Github (thomasp85/patchwork@79223d3)
 pillar             1.6.5       2022-01-25 [1] CRAN (R 4.1.2)
 pkgbuild           1.3.1       2021-12-20 [1] CRAN (R 4.1.2)
 pkgconfig          2.0.3       2019-09-22 [1] CRAN (R 4.1.1)
 pkgload            1.2.4       2021-11-30 [1] CRAN (R 4.1.2)
 png                0.1-7       2013-12-03 [1] CRAN (R 4.1.1)
 prettyunits        1.1.1.9000  2022-01-27 [1] Github (r-lib/prettyunits@8706d89)
 processx           3.5.2       2021-04-30 [1] CRAN (R 4.1.1)
 ps                 1.6.0       2021-02-28 [1] CRAN (R 4.1.1)
 purrr            * 0.3.4.9000  2022-01-27 [1] Github (tidyverse/purrr@5aca9df)
 R.cache            0.15.0      2021-04-30 [1] CRAN (R 4.1.1)
 R.methodsS3        1.8.1       2020-08-26 [1] CRAN (R 4.1.1)
 R.oo               1.24.0      2020-08-26 [1] CRAN (R 4.1.1)
 R.utils            2.11.0      2021-09-26 [1] CRAN (R 4.1.1)
 R6                 2.5.1.9000  2022-01-27 [1] Github (r-lib/R6@1b05b89)
 ragg               1.2.1.9000  2022-01-27 [1] Github (r-lib/ragg@c68c666)
 Rcpp               1.0.8       2022-01-13 [1] CRAN (R 4.1.2)
 RCurl              1.98-1.5    2021-09-17 [1] CRAN (R 4.1.1)
 readr            * 2.1.1.9000  2022-01-30 [1] Github (tidyverse/readr@9610228)
 readxl             1.3.1.9000  2022-01-27 [1] Github (tidyverse/readxl@2ccb82c)
 remotes            2.4.2       2022-01-27 [1] Github (r-lib/remotes@9355549)
 reprex             2.0.1       2021-08-05 [1] CRAN (R 4.1.1)
 rlang              1.0.0.9000  2022-01-30 [1] Github (r-lib/rlang@b0f358e)
 rmarkdown          2.11.12     2022-01-27 [1] Github (rstudio/rmarkdown@b53a7ce)
 rprojroot          2.0.2       2020-11-15 [1] CRAN (R 4.1.1)
 RSQLite            2.2.9       2021-12-06 [1] CRAN (R 4.1.2)
 rstudioapi         0.13.0-9000 2022-01-27 [1] Github (rstudio/rstudioapi@5d0f087)
 Rttf2pt1           1.3.9       2021-07-22 [1] CRAN (R 4.1.1)
 rvest              1.0.2       2021-10-16 [1] CRAN (R 4.1.1)
 S4Vectors        * 0.32.3      2021-11-21 [1] Bioconductor
 sass               0.4.0       2021-05-12 [1] CRAN (R 4.1.1)
 scales             1.1.1.9000  2022-01-29 [1] Github (r-lib/scales@13b01f0)
 sessioninfo        1.2.2       2021-12-06 [1] CRAN (R 4.1.2)
 SparseM          * 1.81        2021-02-18 [1] CRAN (R 4.1.2)
 stringi            1.7.6       2021-11-29 [1] CRAN (R 4.1.2)
 stringr          * 1.4.0.9000  2022-01-27 [1] Github (tidyverse/stringr@85f6140)
 styler           * 1.6.2.9000  2022-01-27 [1] Github (r-lib/styler@9274aed)
 systemfonts        1.0.3       2021-10-13 [1] CRAN (R 4.1.2)
 testthat           3.1.2.9000  2022-01-27 [1] Github (r-lib/testthat@54b9db2)
 textshaping        0.3.6       2021-10-13 [1] CRAN (R 4.1.1)
 tibble           * 3.1.6.9001  2022-01-30 [1] Github (tidyverse/tibble@5d27779)
 tidyr            * 1.1.4.9000  2022-01-27 [1] Github (tidyverse/tidyr@d46d22a)
 tidyselect         1.1.1       2021-04-30 [1] CRAN (R 4.1.1)
 tidyverse        * 1.3.1.9000  2022-01-27 [1] Github (tidyverse/tidyverse@6186fbf)
 topGO            * 2.46.0      2021-10-26 [1] Bioconductor
 tzdb               0.2.0       2021-10-27 [1] CRAN (R 4.1.1)
 usethis            2.1.5.9000  2022-01-27 [1] Github (r-lib/usethis@57b109a)
 utf8               1.2.2       2021-07-24 [1] CRAN (R 4.1.1)
 vctrs              0.3.8       2021-04-29 [1] CRAN (R 4.1.1)
 viridisLite        0.4.0       2021-04-13 [1] CRAN (R 4.1.1)
 vroom              1.5.7.9000  2022-01-30 [1] Github (r-lib/vroom@d91f2aa)
 withr              2.4.3       2021-11-30 [1] CRAN (R 4.1.2)
 xfun               0.29        2021-12-14 [1] CRAN (R 4.1.2)
 xml2               1.3.3       2021-11-30 [1] CRAN (R 4.1.2)
 XVector            0.34.0      2021-10-26 [1] Bioconductor
 yaml               2.2.2       2022-01-25 [1] CRAN (R 4.1.2)
 zlibbioc           1.40.0      2021-10-26 [1] Bioconductor

 [1] /opt/homebrew/lib/R/4.1/site-library
 [2] /opt/homebrew/Cellar/r/4.1.2/lib/R/library

──────────────────────────────────────────────────────────────────────────────